/*
 * SPDX-License-Identifier: BSD-3-Clause
 *
 * Copyright © 2025 Keith Packard
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above
 *    copyright notice, this list of conditions and the following
 *    disclaimer in the documentation and/or other materials provided
 *    with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
 * OF THE POSSIBILITY OF SUCH DAMAGE.
 */

load "test-float.5c"
import Complex;

typedef union {
	complex		c;
	void		i;
	void		n;
} complex_result_t;

complex_result_t complex_nan = { .n = ◊ };

complex_result_t complex_inf = { .i = ◊ };

complex_result_t complex_result(complex z)
{
	return (complex_result_t) { .c = z };
}

complex
complex32(complex x)
{
	return imprecise_arg(x, 24);
}

complex
complex64(complex x)
{
	return imprecise_arg(x, 53);
}

complex
complex80(complex x)
{
	return imprecise_arg(x, 64);
}

complex
complex128(complex x)
{
	return imprecise_arg(x, 113);
}

complex [](complex) complex_convs = { complex32, complex64, complex80, complex128 };

void
print_one_complex(complex z, int i)
{
	print_one_real(creal(z), i);
	printf(", ");
	print_one_real(cimag(z), i);
}

void
print_one_complex_result(complex_result_t z, int i)
{
	union switch (z) {
	case c c:
		print_one_complex(c, i);
		break;
	case i:
		printf("INFINITY, INFINITY");
		break;
	case n:
		printf("NAN, NAN");
		break;
	}
}
void
print_complex(complex[4] x, complex_result_t[4] y)
{
	printf("{ .x = COMPLEX(");
	for (int i = 0; i < 4; i++) {
		if (i > 0)
			printf(", ");
		print_one_complex(x[i], i);
	}
	printf("), .y = COMPLEX(");
	for (int i = 0; i < 4; i++) {
		if (i > 0)
			printf(", ");
		print_one_complex_result(y[i], i);
	}
	printf(") },\n");
}

void
compute_complex_one(complex x, complex(complex) f)
{
	complex[4]		xp;
	complex_result_t[4]	yp;
	for (int i = 0; i < 4; i++) {
		xp[i] = complex_convs[i](x);
		try {
			try {
				yp[i] = complex_result(complex_convs[i](f(imprecise(xp[i], prec))));
			} catch divide_by_zero(real x, real y) {
				if (x == 0)
					raise nan();
				raise infinity(x);
			} catch invalid_argument(string error, int i, poly v) {
				raise nan();
			}
		} catch infinity(real v) {
			yp[i] = complex_inf;
		} catch nan() {
			yp[i] = complex_nan;
		}
	}
	print_complex(xp, yp);
}

void
print_complex_complex(complex[4] x1, complex[4] x2, complex_result_t[4] y)
{
	printf("{ .x1 = COMPLEX(");
	for (int i = 0; i < 4; i++) {
		if (i > 0)
			printf(", ");
		print_one_complex(x1[i], i);
	}
	printf("), .x2 = COMPLEX(");
	for (int i = 0; i < 4; i++) {
		if (i > 0)
			printf(", ");
		print_one_complex(x2[i], i);
	}
	printf("), .y = COMPLEX(");
	for (int i = 0; i < 4; i++) {
		if (i > 0)
			printf(", ");
		print_one_complex_result(y[i], i);
	}
	printf(") },\n");
}

void
compute_complex_complex(complex x1, complex x2, complex(complex, complex) f)
{
	complex[4]		x1p;
	complex[4]		x2p;
	complex_result_t[4]	yp;
	for (int i = 0; i < 4; i++) {
		x1p[i] = complex_convs[i](x1);
		x2p[i] = complex_convs[i](x2);
		try {
			try {
				yp[i] = complex_result(complex_convs[i](f(imprecise(x1p[i], prec),
									  imprecise(x2p[i], prec))));
			} catch divide_by_zero(real x, real y) {
				if (x == 0)
					raise nan();
				raise infinity(x);
			} catch invalid_argument(string error, int i, poly v) {
				raise nan();
			}
		} catch infinity(real v) {
			yp[i] = complex_inf;
		} catch nan() {
			yp[i] = complex_nan;
		}
	}
	print_complex_complex(x1p, x2p, yp);
}
